%% Final Figures for Derivative paper

%% (1) 
clear; % clear all variables
global alpha rho L K I var kappa sig0;
alpha=10; I=10; K=6; sig0=1/I; D=2; var=1; kappa=(1+(I-2)*sig0)/(1-sig0);
rhovec = -1/(K-1)+0.001:0.001:0;
Dvec = [1,2,3,6];
PIPU=zeros(length(rhovec),length(Dvec));PIPUoff=PIPU;errfig1=PIPU;PIa=PIPU;PId=PIPU;
for indD=1:length(Dvec)
    D=Dvec(indD); L=K/D;
    for indrho=1:length(rhovec)
        rho=rhovec(indrho);
        options = optimset('TolX', 10^(-12), 'TolFun', 10^(-12),'MaxFunEvals',10^6,'Display','none');
        [y,fvalue] = fsolve(@Lemma2,alpha/(I-2),options);
        PIPU(indrho,indD)=y;
        PIPUoff(indrho,indD)=y-alpha*var*(1-rho)/(I-2);
        errfig1(indrho,indD)=abs(fvalue);
        PIa(indrho,indD) = alpha*var*(1-rho)/(I-2);
        PId(indrho,indD) = -(PIa(indrho,indD)^2/L/(PIPU(indrho,indD)-alpha*var*(1-rho)/(I-2))+PIa(indrho,indD))/L;
    end
end
    
figure(1)
subplot(1,2,1); plot(rhovec,PIPU(:,2)-alpha*var/(I-2),'-black',rhovec,PIPU(:,3)-alpha*var/(I-2),'-blue',rhovec,PIPU(:,4)-alpha*var/(I-2),'-red','LineWidth',1.2);
xlabel('\rho','FontSize',18); ylabel('$$\hat{\lambda}_k -\lambda_k^c$$','Interpreter','Latex','FontSize',18);
subplot(1,2,2); plot(rhovec,PIPUoff(:,2)-alpha.*var.*rhovec'./(I-2),'-black',rhovec,PIPUoff(:,3)-alpha.*var.*rhovec'./(I-2),'-blue',rhovec,PIPUoff(:,4)-alpha.*var.*rhovec'./(I-2),'-red','LineWidth',1.2);
xlabel('\rho','FontSize',18); ylabel('$$\hat{\lambda}_{k\ell} -\lambda_{k\ell}^c$$','Interpreter','Latex','FontSize',18);



%% (2)
clear('global');
global w1 w2;
alpha=10; I=10; K=4; sig0=1/I; L=2; var=1; kappa=(1+(I-2)*sig0)/(1-sig0);
rho = -0.25; Lvec = [2,4]; wvec = 0:pi/50:pi/2;
Sigma = var.*((1-rho).*eye(K,K)+rho.*ones(K,K));

Lamini = alpha.*diag(diag(Sigma))./(I-2)+eye(K,K);
diff = 1000; count = 0;
while (diff > 10^(-5)) && (count < 30)
	CMat = pinv(Lamini')./(I-1);
	BMat = pinv((1-sig0).*(alpha.*Sigma+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*Sigma);
	Ret1 = (alpha.*Sigma+Lamini)*(BMat*BMat');
	Ret2 = (BMat*BMat');
	Ret1d = diag(diag(Ret1));
	Ret2d = diag(diag(Ret2));
	Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

	diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
	Lamini = Lamnew;
	count = count + 1;
end
diff1 = diff;

if (diff > 10^(-5))
 	Cini = pinv(Lamini')./(I-1);
   	diff = 1000; count = 0;
 	while (diff > 10^(-5)) && (count < 30)
        Laminter = pinv(Cini')./(I-1);
        BMat = pinv((1-sig0).*(alpha.*Sigma+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*Sigma);
      	Ret1 = (alpha.*Sigma+Laminter)*(BMat*BMat');
      	Ret2 = (BMat*BMat');
    	Ret1d = diag(diag(Ret1));
    	Ret2d = diag(diag(Ret2));
        Cnew = pinv(Ret1d*pinv(Ret2d));

       	diff = sqrt(sum(sum((Cini-Cnew).^2)));
     	Cini = Cnew;
    	count = count + 1;
    end    
	diff2 = diff;
else
	diff2 = 10^6;
end    

if (diff1 > diff2) 
  	LamMat = pinv(Cini')./(I-1); 
	err2 = diff2;
else 
	LamMat = Lamini;  
  	err2 = diff1;
end 
    errfig2U= err2; 
if (err2 < 10^(-5)) 
    PIPUU = LamMat(1,1);  
end
        
PIPUk=zeros(length(wvec),length(Lvec));PIPUm=PIPUk;errfig2=PIPUk;
for indL=1:length(Lvec)
    L=Lvec(indL); 
    for indw=1:length(wvec)
        w1 = cos(wvec(indw)); w2=sin(wvec(indw));
        weight=[w1.*ones(1,L),w2.*ones(1,K-L)];
        vard = weight*Sigma*weight';
        weight=[w1.*ones(1,L),w2.*ones(1,K-L)]./sqrt(vard);
    
        if (vard > 10^(-8)) 
             
            WMat=[eye(K,K);weight]; Sigmapl = WMat*Sigma*WMat';   
            Lamini = alpha.*diag(diag(Sigmapl))./(I-2)+eye(K+1,K+1);
            diff = 1000; count = 0;
            while (diff > 10^(-5)) && (count < 30)
                CMat = pinv(Lamini')./(I-1);
                BMat = pinv((1-sig0).*(alpha.*Sigmapl+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*(WMat*Sigma));
                Ret1 = (alpha.*Sigmapl+Lamini)*(BMat*BMat');
                Ret2 = (BMat*BMat');
                Ret1d = diag(diag(Ret1));
                Ret2d = diag(diag(Ret2));
                Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

                diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
                Lamini = Lamnew;
                count = count + 1;
            end
            diff1 = diff;

            if (diff > 10^(-5))
                Cini = pinv(Lamini')./(I-1);
                diff = 1000; count = 0;
                while (diff > 10^(-5)) && (count < 30)
                    Laminter = pinv(Cini')./(I-1);
                    BMat = pinv((1-sig0).*(alpha.*Sigmapl+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*(WMat*Sigma));
                    Ret1 = (alpha.*Sigmapl+Laminter)*(BMat*BMat');
                    Ret2 = (BMat*BMat');
                    Ret1d = diag(diag(Ret1));
                    Ret2d = diag(diag(Ret2));
                    Cnew = pinv(Ret1d*pinv(Ret2d));

                    diff = sqrt(sum(sum((Cini-Cnew).^2)));
                    Cini = Cnew;
                    count = count + 1;
                end    
                diff2 = diff;
            else
                diff2 = 10^6;
            end    

            if (diff1 > diff2) 
                LamMat = pinv(Cini')./(I-1); 
                err2 = diff2;
            else 
                LamMat = Lamini;  
                err2 = diff1;
            end 
            errfig2(indw,indL)= err2+errfig2U; 
            if (err2 < 10^(-5))
                PIPUMat = inv([eye(K,K);weight]'*inv(LamMat)*[eye(K,K);weight]);
                PIPUk(indw,indL) = PIPUMat(1,1);
                PIPUm(indw,indL) = PIPUMat(min(K,L+1),min(K,L+1));  
                
            end 
        end
    end      
end

figure(2)
subplot(1,2,1); plot(cos(wvec),PIPUk(:,1)-alpha*var/(I-2),'-black',cos(wvec),zeros(1,length(wvec)),'--black',cos(wvec),(PIPUU-alpha*var/(I-2)).*ones(length(wvec)),':black','LineWidth',1.2);
xlabel('w_k','FontSize',18); ylabel('$$\hat{\lambda}_k-\lambda^c_k$$','Interpreter','Latex','FontSize',18);
xlim([min(cos(wvec)),max(cos(wvec))]);
colormap('gray');
subplot(1,2,2); plot(cos(wvec),PIPUm(:,1)-alpha*var/(I-2),'-black',cos(wvec),zeros(1,length(wvec)),'--black',cos(wvec),(PIPUU-alpha*var/(I-2)).*ones(length(wvec)),':black','LineWidth',1.2);
xlabel('w_k','FontSize',18); ylabel('$$\hat{\lambda}_{m}-\lambda^c_m$$','Interpreter','Latex','FontSize',18);
xlim([min(cos(wvec)),max(cos(wvec))]);
colormap('gray');
 

%% (3)
clear;clc;clear('global');
global alpha Sigma0 I qvec sig0; 
alpha=10; I=10; K=3; sig0=1/I; var=1; kappa=(1+(I-2)*sig0)/(1-sig0); L=1;
rhol=-0.1; ql=10; 
rhohvec = [-0.201:0.01:-0.01,0.01:0.01:0.201]; qhvec=[-2:0.1:2].*ql;  
errcount = 0;
for indrho = 1:length(rhohvec)
    Sigma0 = var.*[1,rhohvec(indrho),rhohvec(indrho);
                   rhohvec(indrho),1,rhol;
                   rhohvec(indrho),rhol,1]; 
    % uncontingent
    Lamini = alpha.*diag(diag(Sigma0))./(I-2)+eye(K,K);
	diff = 1000; count = 0;
    while (diff > 10^(-5)) && (count < 30)
    	CMat = pinv(Lamini')./(I-1);
     	BMat = pinv((1-sig0).*(alpha.*Sigma0+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*Sigma0);
    	Ret1 = (alpha.*Sigma0+Lamini)*(BMat*BMat');
     	Ret2 = (BMat*BMat');
     	Ret1d = diag(diag(Ret1));
       	Ret2d = diag(diag(Ret2));
     	Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

      	diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
    	Lamini = Lamnew;
      	count = count + 1;
    end
  	diff1 = diff;

   	if (diff > 10^(-5))
    	Cini = pinv(Lamini')./(I-1);
      	diff = 1000; count = 0;
      	while (diff > 10^(-5)) && (count < 30)
         	Laminter = pinv(Cini')./(I-1);
         	BMat = pinv((1-sig0).*(alpha.*Sigma0+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*Sigma0);
        	Ret1 = (alpha.*Sigma0+Laminter)*(BMat*BMat');
           	Ret2 = (BMat*BMat');
         	Ret1d = diag(diag(Ret1));
           	Ret2d = diag(diag(Ret2));
          	Cnew = pinv(Ret1d*pinv(Ret2d));

         	diff = sqrt(sum(sum((Cini-Cnew).^2)));
         	Cini = Cnew;
          	count = count + 1;
        end    
      	diff2 = diff;
    else
     	diff2 = 10^6;
    end    

  	if (diff1 > diff2) 
     	LamMatU = pinv(Cini')./(I-1); 
    	errU = diff2;
    else 
     	LamMatU = Lamini;  
        errU = diff1;
    end 
        
    bestDini = -sign(rhohvec(indrho));
    for indq = 1:length(qhvec)
        qvec = [qhvec(indq),ql,ql];   
        utilC(indrho,indq)=qvec*(alpha.*Sigma0)*qvec'*(0.5*I*(I-2)/(I-1)/(I-1));
        
        if (errU < 10^(-5))
            utilU(indrho,indq) = qvec*((alpha.*Sigma0)*pinv(alpha.*Sigma0+LamMatU')*(0.5.*alpha.*Sigma0+0.5.*LamMatU+0.5.*LamMatU')*pinv(alpha.*Sigma0+LamMatU)*(alpha.*Sigma0))*qvec';
            maxutil(indrho,indq)=max(utilC(indrho,indq),utilU(indrho,indq));
        else
            utilU(indrho,indq)=utilC(indrho,indq);
            maxutil(indrho,indq)=max(utilC(indrho,indq),utilU(indrho,indq));
        end
        
        % best derivative
        options = optimset('TolX', 10^(-12), 'TolFun', 10^(-12),'MaxFunEvals',10^6,'Display','none');
        [bestD,fvalue] = fminunc(@Ex4,bestDini,options);
        
        weight = [bestD,1,1];
        vard=weight*Sigma0*weight'; weight = weight./sqrt(vard);
        WMat = [eye(K,K);weight]; Sigmapl = WMat*Sigma0*WMat';
        Lamini = alpha.*[diag(diag(Sigma0)),zeros(K,L);zeros(L,K),alpha*Sigma0(1,1).*eye(L,L)]./(I-2)+eye(K+L,K+L);
        diff = 1000; count = 0;
        while (diff > 10^(-5)) && (count < 30)
            CMat = pinv(Lamini')./(I-1);
            BMat = pinv((1-sig0).*(alpha.*Sigmapl+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*(WMat*Sigma0));
            Ret1 = (alpha.*Sigmapl+Lamini)*(BMat*BMat');
            Ret2 = (BMat*BMat');
            Ret1d = diag(diag(Ret1));
            Ret2d = diag(diag(Ret2));
            Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

            diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
            Lamini = Lamnew;
            count = count + 1;
        end
        diff1 = diff;

        if (diff > 10^(-5))
            Cini = pinv(Lamini')./(I-1);
            diff = 1000; count = 0;
            while (diff > 10^(-5)) && (count < 30)
                Laminter = pinv(Cini')./(I-1);
                BMat = pinv((1-sig0).*(alpha.*Sigmapl+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*(WMat*Sigma0));
                Ret1 = (alpha.*Sigmapl+Laminter)*(BMat*BMat');
                Ret2 = (BMat*BMat');
                Ret1d = diag(diag(Ret1));
                Ret2d = diag(diag(Ret2));
                Cnew = pinv(Ret1d*pinv(Ret2d));

                diff = sqrt(sum(sum((Cini-Cnew).^2)));
                Cini = Cnew;
                count = count + 1;
            end    
            diff2 = diff;
        else
            diff2 = 10^6;
        end    

        if (diff1 > diff2) 
            LamMat = pinv(Cini')./(I-1); 
            err2 = diff2;
        else 
            LamMat = Lamini;  
            err2 = diff1;
        end 
        if (err2 < 10^(-5))
            Sighalf = WMat*Sigma0; 
            utilDeriv(indrho,indq) = qvec*(alpha.*Sigma0)*WMat'*pinv(alpha.*Sigmapl+LamMat')*(0.5.*alpha.*Sigmapl+0.5.*LamMat+0.5.*LamMat')*pinv(alpha.*Sigmapl+LamMat)*WMat*(alpha.*Sigma0)*qvec';
            wratio(indrho,indq) = weight(1)/weight(2);
            errfig4(indrho,indq) = err2;
        else
            utilDeriv(indrho,indq) = max(utilC(indrho,indq),utilU(indrho,indq))-0.25;
            wratio(indrho,indq) = 0;
            errfig4(indrho,indq) = err2;
        end
    end
end

figure(4)
surf(qhvec./ql,rhohvec./rhol,utilDeriv-maxutil);
xlabel('q_k^+/q_l^+','FontSize',14);
ylabel('\rho_k^+/\rho_l^+','FontSize',14);
zlabel('Welfare Change','FontSize',14);
colormap('gray');

figure(42)
surf(qhvec./ql,rhohvec./rhol,wratio);
xlabel('q_k^+/q_l^+','FontSize',14);
ylabel('\rho_k^+/\rho_l^+','FontSize',14);
zlabel('\omega_k/\omega_l','FontSize',14);
colormap('gray');

%% (4)
clear;clc; clear('global');
global rho L K rhopl;
alpha=10; I=10; K=3; L = 2; sig0=1/I; var=1; kappa=(1+(I-2)*sig0)/(1-sig0); Z = K+L;
rhovec = [-0.2:0.01:0.2];
qvec = [10.*ones(1,K),1.*ones(1,L)];

utilDiff = zeros(length(rhovec),length(rhovec));

for indrho = 1:length(rhovec)
    rho=rhovec(indrho);
    Sigma0 = (1-rho).*eye(Z,Z)+rho.*ones(Z,Z);
    
    WMat = eye(K,Z); Sigmapl = WMat*Sigma0*WMat';
  	Lamini = alpha.*diag(diag(Sigmapl))./(I-2)+eye(K,K);
    diff = 1000; count = 0;
    while (diff > 10^(-5)) && (count < 30)
        CMat = pinv(Lamini')./(I-1);
     	BMat = pinv((1-sig0).*(alpha.*Sigmapl+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*(WMat*Sigma0));
      	Ret1 = (alpha.*Sigmapl+Lamini)*(BMat*BMat');
        Ret2 = (BMat*BMat');
     	Ret1d = diag(diag(Ret1));
      	Ret2d = diag(diag(Ret2));
     	Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

        diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
    	Lamini = Lamnew;
     	count = count + 1;
    end
	diff1 = diff;

	if (diff > 10^(-5))
     	Cini = pinv(Lamini')./(I-1);
       	diff = 1000; count = 0;
     	while (diff > 10^(-5)) && (count < 30)
         	Laminter = pinv(Cini')./(I-1);
          	BMat = pinv((1-sig0).*(alpha.*Sigmapl+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*(WMat*Sigma0));
          	Ret1 = (alpha.*Sigmapl+Laminter)*(BMat*BMat');
          	Ret2 = (BMat*BMat');
         	Ret1d = diag(diag(Ret1));
          	Ret2d = diag(diag(Ret2));
           	Cnew = pinv(Ret1d*pinv(Ret2d));

        	diff = sqrt(sum(sum((Cini-Cnew).^2)));
        	Cini = Cnew;
          	count = count + 1;
        end    
        diff2 = diff;
    else
        diff2 = 10^6;
    end    

 	if (diff1 > diff2) 
     	LamMat = pinv(Cini')./(I-1); 
      	err2 = diff2;
    else 
     	LamMat = Lamini;  
      	err2 = diff1;
    end
    if (err2 < 10^(-5))
        utilU(indrho) = qvec*((alpha.*Sigma0)*WMat'*pinv(alpha.*Sigmapl+LamMat)*(0.5.*alpha.*Sigmapl+0.5.*LamMat+0.5.*LamMat')*pinv(alpha.*Sigmapl+LamMat)*WMat*(alpha.*Sigma0))*qvec';
     	errfig5U(indrho) = err2;
    end
            
            
    for indrhopl = 1:length(rhovec)
        rhopl=rhovec(indrhopl)*1.2; err2 = 0;
        Sigma0 = (1-rho).*eye(Z,Z)+rho.*ones(Z,Z);

        options = optimset('TolX', 10^(-12), 'TolFun', 10^(-12),'MaxFunEvals',10^6,'Display','none');
        [y,fvalue] = fsolve(@FWC,0,options);
        c = y; b = 1/sqrt(1-rho); a = (rhopl-(b+L*c)*rho)/(1+(K-1)*rho);
        WMat = [eye(K,Z);a.*ones(L,K), b.*eye(L,L)+c.*ones(L,L)];
        Sigmapl = WMat*Sigma0*WMat';
        
    	if (abs(fvalue) < 10^(-5))
            Lamini = alpha.*diag(diag(Sigma0))./(I-2)+eye(Z,Z);
            diff = 1000; count = 0;
            while (diff > 10^(-5)) && (count < 30)
                CMat = pinv(Lamini')./(I-1);
                BMat = pinv((1-sig0).*(alpha.*Sigmapl+Lamini)+sig0.*(I-1).*Lamini')*(alpha.*(WMat*Sigma0));
                Ret1 = (alpha.*Sigmapl+Lamini)*(BMat*BMat');
                Ret2 = (BMat*BMat');
                Ret1d = diag(diag(Ret1));
                Ret2d = diag(diag(Ret2));
                Lamnew = (Ret1d*pinv(Ret2d))'./(I-1);

                diff = sqrt(sum(sum((Lamini-Lamnew).^2)));
                Lamini = Lamnew;
                count = count + 1;
            end
            diff1 = diff;

            if (diff > 10^(-5))
                Cini = pinv(Lamini')./(I-1);
                diff = 1000; count = 0;
                while (diff > 10^(-5)) && (count < 30)
                    Laminter = pinv(Cini')./(I-1);
                    BMat = pinv((1-sig0).*(alpha.*Sigmapl+Laminter)+sig0.*(I-1).*Laminter')*(alpha.*(WMat*Sigma0));
                    Ret1 = (alpha.*Sigmapl+Laminter)*(BMat*BMat');
                    Ret2 = (BMat*BMat');
                    Ret1d = diag(diag(Ret1));
                    Ret2d = diag(diag(Ret2));
                    Cnew = pinv(Ret1d*pinv(Ret2d));

                    diff = sqrt(sum(sum((Cini-Cnew).^2)));
                    Cini = Cnew;
                    count = count + 1;
                end    
                diff2 = diff;
            else
                diff2 = 10^6;
            end    

            if (diff1 > diff2) 
                LamMat = pinv(Cini')./(I-1); 
                err2 = diff2;
            else 
                LamMat = Lamini;  
                err2 = diff1;
            end 
            if (err2 < 10^(-5))
                Lam(:,:,indrho,indrhopl)=LamMat;
                utilZ(indrho,indrhopl) = qvec*((alpha.*Sigma0)*WMat'*pinv(alpha.*Sigmapl+LamMat)*(0.5.*alpha.*Sigmapl+0.5.*LamMat+0.5.*LamMat')*pinv(alpha.*Sigmapl+LamMat)*WMat*(alpha.*Sigma0))*qvec';
                errfig5(indrho,indrhopl) = err2;
                utilDiff(indrho,indrhopl)=utilZ(indrho,indrhopl)-utilU(indrho);
            end
        end
    end
end

figure(5) 
surf(rhovec,rhovec.*1.2,utilDiff);
ylabel('\rho','FontSize',14); 
xlabel('\rho^+','FontSize',14);
zlabel('Welfare Change','FontSize',14);
colormap('gray');


